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ABSTRACT 

We attempt to recover the mean vertical velocity and vertical velocity dispersion as a 
function of the Galactic height for a sample draw n from a realistic Galaxy d istribution 
function by following the method presented in Moni Bidi n et all (l2Q12al) . We find 
that, for the sample size used, the observational error in the vel ocities is much sm aller 
than the Poisson noise which has not been accounted for by Mo ni Bidin et all We 
repeat the analysis on a large number of samples to estimate the contribution of 
the Poisson noise and to uncover any systematics. We find that the dispersion is 
systematically overestimated at low Galactic heights and slightly underestimated at 
high Galactic heights leading to an underestimate of the gradient of the dispersion 
with Galactic height. The causes of t he systemati c s are revealed by repeating the 
calculation using a method inspired bv lGirard et al.l (|2QQ6f ). This method recovers the 
expected dispersion much more successfully and in particular yields a gradient of the 
dispersion with Galactic height which is ap proximately three times that found using 
the method presented bv lMoni Bidin et al.l 

Key words: methods: data analysis - methods: numerical - Galaxy: kinematics and 
dynamics - Galaxy: solar neighbourhood 



1 INTRODUCTION 



iMoni Bidin et all (|2012bh used the Jeans equation to con- 
strain the stellar mass density at heights 1.5 kpc ^ z ^ 4kpc 
from the Galactic plane using a sample of 412 red giants 
and concluded that there is a lack of dark ma tter in the so- 
lar neighbourhood. iBovv &; Tremainel (|2012l) have pointed 
out tha t one of the assumptions made bv lMoni Bidin et al.l 
(|2012bh in their Jeans-equation analysis (specifically the as- 
sumption that the mean azimuthal velocity is independent 
of Galactocentric radius at all heights) is false and a reanal- 
ysis of the data without this assumption leads to a non-zero 
local dark matter density. Despite the fact that the data 
now appear to conform with standard expectation, it is wise 
to check that all parts of the analysis are sound before the 
matter is put to rest. To use the Jeans equation with con- 
fidence one needs to be sure that an underlying population 
has been isolated and that the mean velocities and veloc- 
ity dispersions of that population can be r eliably calculated. 
Here w e investigate the method used by IMoni Bidin et al.l 
(2012a|, hereafter MB) to calculate the mean velocities and 
velocity dispersions of the thick disc. It is these qu anti- 
ties w hich ar e then used by IMoni Bidin et al ] (|2012bh and 
IBovv fe Tremainel (2012) in their Jeans-equation analysis, 
so it is crucial that they are calculated correctly and that 
their associated errors are realistic. 
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MB use the probability plot method which we detail in 
Section [2] We proceed by producing a sampl e drawn from 
the realistic Galaxy distribution functions of lBinnevI (|2012ah 
chosen to be similar to the sample of MB. These distribu- 
tion functions are briefly discussed in Section [3] In Section [H 
we attempt to recover the known mean velocity and veloc- 
ity dispersions of the thick disc using the method presented 
by MB. We also implement the error analysis used by MB, 
which does not include any estimates of the Poisson noise. 
In Section 14.1.21 the analysis is repeated for a large number 
of samples to investigate the effects of the Poisson noise and 
to uncover any systematics introduced by the method. Fi- 
nally in Section [5l t he same data s et is a nalysed by a similar 
method inspired bv lGirard et al.l (|2006h and the results are 
compared. 



2 PROBABILITY PLOT METHOD 

The data analysis used by MB uses the probability plot to 
determine the mean and standard deviation of a sample. 
Here we briefly present the method and give a simple exam- 
ple to demonstrate its use. 

Suppose we have N ordered data points. If these N data 
points are drawn from a normal distribution of mean fi and 
standard deviation a then, for large N, the ith data point, 
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Xi, should satisfy 



F(xi]fi,a) 



N + l' 



F(x;{i, a) is the cumulative distribution function for 
mal distribution, which is given by 

F( W ) = I[l +erf (^)], 



(1) 



(2) 



where erf is the error function. Hence the ith data point 
should lie at Ci standard deviations from the mean where 



(3) 



Therefore a plot of x% against d should have gradient a and 
intercept /x, which may be found by linear regression. Such 
a plot is termed a probability plot. 

Clearly this approach only works exactly if N is large 
and the data have been drawn from a single normal distri- 
bution. However we can use it to estimate the mean and 
standard deviation of a sample drawn from any underly- 
ing distribution which is approximately Gaussian. Here we 
show how the method operates for a simple case. We draw 50 
data points from a normal distribution of mean fi — and 
standard deviation a — 1, and assess how well the above 
method can recover these quantities. Fig. [T] shows the result 
for one randomly drawn sample. The measured standard 
deviation for this sample is s — 1.006. From the linear re- 
gression the mean is estimated asm = 0.156 ±0.141 and the 
standard deviation s — 1.075 =b 0.152 where the errors are 
given by the deviations of the points away from a straight 
line. With N data points the expected error in the mean is 
^/cr 2 JN and the expected error in the variance is approx- 
imately y^2a 4 /N. Therefore, for this sample we expect an 
error in the mean of 0.141 and an error in the standard de- 
viation of 0.1. This is well represented by the errors in the 
linear regression. Thus the method does work well when the 
data are drawn from an underlying Gaussian distribution 
and the errors from the linear regression are comparable to 
the expected Poisson noise. 

When the data have been drawn from the sum of two 
Gaussian distributions, we can estimate the means and stan- 
dard deviations by fitting straight line segments to differ- 
ent parts of the plot. This is the procedure followed by 
iBochanski et al.l ([2007) to calculate the velocity dispersions 
of the thin and thick discs using a tracer population of 7398 
M dwarfs. 



3 DYNAMICAL GALAXY MODELS 

To test the method used by MB we require a sample drawn 
from a realistic Galaxy distribution function for which we 
kn ow the u nderly ing velocity moments. We use the models 
of iBinnevI (l2012ah . Thes e models were developed and dis - 
cussed by Binney (|2010h and iBinnev fe McMillan! (|201lh . 
The distribution function is a function of the actions in an 
axisymmetric potential: the radial action J r , the vertical ac- 
tion J z and the ^-component of the angular momentum, L z . 
The advantage of this approach is that the distribution func- 
tion clearly satisfies the Jeans' theorem as the actions are 
isolating integrals. These models consist of a thick and thin 
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Figure 1. Demonstration of finding the mean and standard devi- 
ation of a sample. The top panel shows a histogram of data drawn 
from a normal distribution of mean zero and unit standard de- 
viation shown by the line. The bottom panel shows the ordered 
data plotted against the expected deviation in units of the stan- 
dard deviation, a, given by Eq. [3] A straight line fit to this plot 
yields an estimate for the standard deviation as the gradient and 
an estimate for the mean as the intercept. For this example the 
standard deviation is estimated as 1.075 and the mean as 0.156. 



disc composed of quasi-isothermal distribution functions: 



f(J r ,Jz,L z ) — f ar ( J r , L z )f az [J z 



where 
and 



[1 +tanh(L 2 /L )]e" 



tJ r / 'a 



(4) 



(5) 



(6) 



Here R C (L Z ) is the radius of a circular orbit with z- 
component of angular momentum, L z . k(L z ), v(L z ) and 
Q(L Z ) are the radial, vertical and circular epicycle frequen- 
cies respectively and E(L^) is the approximate surface den- 
sity of the disc. The factor of [1 + tanh(L 2 /Lo)] eliminates 
retrograde stars. a r (L z ) and a z (L z ) are exponentially de- 
caying function of R c which control the radial and vertical 
velocity dispersions: 



(Tr(L z ) 

cr z {L z ) 



<r r oe 



q(R -R c )/R d 
q(R -R c )/R d 



(7) 
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Table 1. Parameters for the Binney distribution function used 
throughout the paper. Velocity dispersions have units kms -1 and 
scale lengths have units kpc. 
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where a r o and a z o are approximately equal to the radial and 
vertical velocity dispersions at the solar radius, Ro, and Rd 
is the scale length of the disc. 

The thick disc consists of a single quasi- isothermal 
distribution function, /thick, of fixed age whilst the thin 
disc is a superposition of quasi-isothermal distribution func- 
tions of differing a pes from zero to r m = 10 Gyr. Following 
lAumer &; Binnevl (|2009h we also adopt an age dependence 
for the velocity dispersion such that 



a r (L z ,r] 
ct z (L z ,t] 



( 1 

= (J rO — 



T + Ti 



T + T\ 



( T + Tl \ 
<TzO ( ) 



^q(R -R c )/R d 



a(R -R c )/R d 



(8) 



where we set f3 = 0.33 and n = 0.11 Gyr. We also assume a 
decreasing star formation rate with time with characteristic 
time to = 8 Gyr such that the full distribution for the thin 
disc is given by 



/ti 



(9) 



to(eWto-i) 

The ratio of the thick to thin disc distribution functions is 
controlled by a parameter, k t hk = 0.224, which implies a 
fraction 0.224/1.224 = 0.18 of the disc stars belong to the 
thick disc and implies a Solar neighbourhood ratio of 0.28. 
We set the parameters to the values given in Table [T] The 
density profile and velocity dispersion for this distribution 
function as a function of Galactic height are shown in Fig. [2] 
We limit our investigation to just the W component 
of the velocity (the component along the z direction). We 
are not interested in the full distribution function so we 
marginalise over the other two velocity components (U and 
V) to find the number of stars per unit W velocity per unit 
volume as 



nw{R, 2, W) 



/oo rc 
-oo J —c 



dVf{J r ,J z ,L z 



(10) 



where in practice the limits of the integrals are finite as 
the distribution function falls off rapidly at large velocities. 
The transformation from polar positions and velo cities to 
actions is carried out by the algorithm presented by iBinnevI 
_(l2012bh. We use an adjusted version of Potential II from 
iDehnen &; Binne vl (jl998h which consists of a thin and thick 
disc, a gas disc and two spheroids representing the bulge and 
the halo. We have increased the scale-height of the thin disc 
to 360 pc and increased the mass of the thin disc such that 
the circular velocity at the solar radius is 220 kms -1 . 
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Figure 2. Density profile and velocity dispersion for the Binney 
distribution function used in this paper as a function of Galactic 
height. The dashed line gives the contribution from the thin disc 
and the dotted line from the thic k disc. The data poin t s are taken 
from lGilmore Reidl (119831 ) and lKuiiken Gilmore] (119891 ). 



From this distribution function we are able to draw 
a sample of stars. The sample is selected by following a 
rejection algorithm. If we wish to draw stars which all 
lie at the solar radius, Ro, and which lie in the ranges 
z m in ^ z ^ 2 max and \W\ ^ Wmax then we first note that 
the maximum of the distribution function in this range oc- 
curs at z = Zmin and W = 0. This gives us a normalisa- 
tion. We then proceed by drawing trial values of z = zt 
and W = Wt from uniform distributions over the required 
ranges and accepting this trial as a data point with proba- 
bility n w (Ro,z t ,W t )/ n w (Ro , ^min , 0) . 

The obvious advantage of drawing sample data from 
a known distribution function is that we know exactly the 
underlying properties of the distribution. In this paper we 
focus on calculating the mean W velocity and the W velocity 
dispersion of the thick disc, a z . As the distribution function 
is a symmetric function of W we expect (W) = 0. a z is given 
by 



a 2 z {R,z)= J - 



^ dUJZy dVj^ dWW 2 f t ^(J r , J Z ,L Z ) 
J^dUf^dVj^dWf t ^(J r ,J z ,L z ) ■ 



(11) 



4 MB METHOD 

MB use a sample of 412 red giants which lie in the range 
1.3 kpc < z < 5kpc and \W\ < 150 kms" 1 . The cut in W 
is performed to remove contamination from the halo. MB 
estimate an error of 0.7 kms -1 in the radial velocity mea- 
surements and an error of approximately 20 per cent in the 
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Figure 3. Sample of 412 stars drawn from the Binney distribu- 
tion function that emulates the MB sample. 

distance^]. We use the procedure outlined above to draw 412 
stars from the Binney distribution function which all lie at 
the solar radius and inside the range probed by MB. As the 
W velocity is nearly entirely radial velocity error we include 
a random Gaussian error of 0.7 km s -1 to the W velocities 
and we assume that the full distance error of 20 per cent 
corresponds to a 20 per cent error in the z values. In order 
to correctly account for stars which may have entered our 
sample due to the error in their distances we increase the 
sample range to 1 kpc ^ z ^ 6 kpc and then cut out any 
stars which lie outside 1.3 kpc < z ^ 5 kpc after the error 
has been included. Histograms of the resulting sample are 
shown in Fig. [3] 

We now follow the same procedure as MB to extract 
the mean W velocity and the W velocity dispersion, a z . We 
first bin the data in z with bin centres spaced by 0.1 kpc in 
the range 1.5 kpc ^ z ^ 4.5 kpc. The bin sizes are allowed 
to vary such that we have 100 data per bin for z ^ 2.1 kpc, 
80 data per bin for 2.2 kpc ^ z ^ 2.4 kpc and 50 data per 
bin for z ^ 2.5 kpc. For each binned subset of the sample 
we follow the method outlined in Section [2] by ordering the 
sample in W and performing a linear regression between 
the sample velocities and the expected deviations to find 
the mean velocity of the bin and the velocity dispersion. 

For z ^ 2.5 kpc we expect a non-negligible thin disc 
contamination. MB try to isolate the thick disc contribu- 
tion by only fitting the wings of the distribution where, as 
the thick disc has a higher velocity dispersion than the thin 
disc, the data are assumed to be contributed by thick disc 
stars. Therefore, we first sort all the data in the bin and 
assign each an expected deviation but only fit a straight 



1 MB state that there is an additional 10 — 20 per cent systematic 
error in the distances to thin disc stars due to the thin disc stars 
not following the assumed age and metallicity distributions. We 
ignore this additional error here. This error will increase thin disc 
contamination at low z but should not affect the determination of 
the velocity dispersion or conclusions presented here significantly. 



line to the data which have \W\ > 30kms -1 . Each wing is 
fitted separately. MB do not make it completely clear how 
they combine the fits of each wing but here we adopt the 
procedure of fitting each wing independently and then cal- 
culating the mean and standard deviation by an average of 
the intercepts and gradients respectively. 

MB also ignore any points which seem to be outliers in 
the probability plots when performing the linear regression. 
We simulate this effect by ignoring the most negative and 
most positive data point when fitting a straight line to the 
probability plot. 



4.1 Errors 

4-1.1 Observational Errors 

MB estimate the errors in their calculated W moments by 
essentially only considering the observational error in the 
W velocity as follows. In each bin MB add random Gaus- 
sian errors for the distance and radial velocity to each data 
point to generate 1000 samples but do not re-bin the data 
at all. The errors are then estimated as the standard devia- 
tion of the estimates obtained from each sample. Following 
MB we take the W velocity error to be 0.7 km s -1 . The re- 
sults of this procedure are shown in Fig. 2] The error bars 
are very small giving the impression we have very precise 
results. However the data are clearly scattered around the 
true result by amounts much greater than the error bars. 
This is because we have ignored two much larger sources of 
error: the errors in the distances moving stars from bin to 
bin and the Poisson noise. 

MB estimate the distance error to be approximately 
20 per cent. As well as a W velocity error, we add a random 
Gaussian error of 20 per cent to the z coordinates of the 
data, re-bin the data and recalculate the mean velocities 
and dispersions. Repeating this 1000 times we calculate the 
errors as the standard deviations of the estimates. These 
results are shown in Fig. [5] The observational errors are now 
much larger and the results are consistent with the truth 
within the errors. We have not yet made any estimate for 
the Poisson noise of the estimate but it seems that, as the 
data are consistent with the truth, the observational errors 
are of the same order as the Poisson noise. 

One standard deviation lines are shown in Fig. [5] 
Adding errors in z to the data shifts stars from bin to bin. 
With such a large distance error we have many samples in 
which higher velocity stars have been pulled down to lower 
Galactic heights and lower velocity stars are displaced to 
greater heights. This has the effect of flattening the velocity 
distribution when averaged over the many samples. 



4-1.2 Poisson Noise and systematics 

With such a small sample it is difficult to disentangle the 
Poisson noise from systematic errors arising from the MB 
procedure. However, as we have direct access to the dis- 
tribution function we can estimate the Poisson noise by re- 
peatedly drawing samples from the distribution function and 
evaluating the observables for each sample. Therefore we 
draw 100 samples of 412 stars and repeat the above proce- 
dure for each sample. We then estimate the average sample 
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Figure 4. Mean velocity and W velocity dispersion, a z , 
against height above Galactic plane. The points give the values 
calculated in each bin using the MB method along with the errors 
estimated by 1000 samples adding random errors of O.Tkms- 1 to 
the W velocity. The red dashed lines show the exact mean velocity 
and velocity dispersion calculated from the Binney distribution 
function from which the sample was drawn. Clearly the error bars 
do not give a good estimate of the deviation from the true value. 

mean and dispersion in each bin along with the Poisson er- 
ror in both quantities by calculating the standard deviations. 
The results are shown in Fig. [6] 

The calculation of the mean W velocity is entirely con- 
sistent with being zero as required, but the recovery of the W 
velocity dispersion curve is less successful. For low Galactic 
heights we are overestimating the velocity dispersion whilst 
for larger Galactic height we are slightly underestimating 
the dispersion. 

For low z we are ignoring all data for which \W\ < 
30 km s _1 when fitting a straight line to the data. This means 
we give more weight to data with higher W velocities and 
so the distribution seems broader than it actually is. The 
probability plot is particularly sensitive at the wings. If we 
consider an ordered data set that is drawn from a known 
underlying Gaussian distribution, we can assign each a value 
of Ci by the method outlined in Section [2] If we add a single 
point which is smaller than all the other data points but 
still drawn from the underlying distribution the probability 
of the new point lying above the line with correct mean 
and standard deviation, but still lower than its neighbouring 
data point, is equal to the probability of it lying beneath 
the line. As there is a much larger range of values below 
the line than above, the estimated gradient in this region 
will in general be overestimated. We need to use a sufficient 
number of stars to perform the linear fit in order to reduce 
this effect. 

At high z we have very few stars in the sample so in 
order to fill the bin with enough stars we must include stars 
at lower z. In general these stars have smaller velocities and 



Figure 5. Mean W velocity and W velocity dispersion, a z , 
against height above Galactic plane for the sample of 412 stars. 
The points give the values calculated in each bin using the MB 
method. The dotted lines show the one standard deviation limits. 
These were estimated by calculating the standard deviations of 
the calculated values for 1000 samples formed by adding random 
errors of 0.7 km s -1 to the W velocity and a 20 per cent error on 
the distance with re-binning. The red dashed lines show the exact 
mean velocity and velocity dispersion calculated from the Binney 
distribution function from which the sample was drawn. 

so the resulting velocity dispersion for the bin is reduced. A 
very minor effect may also be due to removing stars which 
have \W\ > 150kms _1 to avoid halo contaminants, so the 
distribution is unfairly weighted by low-velocity stars and 
the dispersion is underestimated. From Fig. [6] we can per- 
form a simple linear fit to the (z,a z ) plot to find that the 
data points imply a gradient of da z /dz = 2.25 kms -1 kpc -1 . 
A linear fit to the true dispersion curve gives a gradient of 
da z /dz = 6.82 kins" 1 kpc -1 so the MB method underesti- 
mates the gradient by a factor of three. 



5 COMPARISON WITH OTHER WORK 

MB state that their gradient of <j z with z is shallower than 
previous authors' work. The s ample studied by M B is a sub- 
set of the sample studied by iGirard et al. (2006, hereafter 
G06). These authors found a gradient of the V velocity dis- 
persion a factor of two higher than MB, which MB claim is 
due to G06 not removing thin disc and halo contaminants in 
their analysis. We conclude by following a method inspired 
by the method used by G06 to test whether their results are 
more secure than those reported by MB. The G06 method is 
very similar to that used by MB. The authors have a sample 
of approximately 1200 stars. They split each data point into 
100 subunits to account for the distance error and form bins 
of 100 subunits to estimate the U and V velocity dispersion. 
The probability plot method is used to estimate the disper- 
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Figure 6. MB method: Average mean W velocity and W ve- 
locity dispersion, cr 2 , against height above Galactic plane for the 
100 samples of 412 stars. The points give the average values calcu- 
lated in each bin using the MB method. The error bars show the 
standard deviation of the calculated values in each bin. The red 
dashed lines show the velocity moments of the thick disc calcu- 
lated using the underlying distribution function. The dispersion 
is overestimated at low z and slightly underestimated at high z 
leading to an underestimated gradient. 

sion but crucially only the central 80 per cent of the data is 
used in the linear fit and the central region is not excluded 
for any of the bins. We follow a similar, but simpler, method 
on each of the 100 samples of 412 stars. We do not split the 
data points into subunits. We use bins spaced by 0.1 kpc 
containing 32 data points, as we have only a third of the 
number of data points in the G06 sample, and we use only 
the central 80 per cent of the data in each bin for the linear 
fit in the probability plot method. Fig. [3 shows the result of 
this experiment. 

The estimate of the velocity dispersion provided by the 
G06 method is more reliable than the MB method. At low z 
the dispersion is now marginally underestimated which is to 
be expected due to the thin disc contamination. However, 
even with fewer stars in the bin, the error in the dispersion 
at low z is smaller than the equivalent errors in the MB 
method. This is a clear reflection of the dangers of only using 
the wings of the distribution to calculate the dispersion. At 
high z the problem of preferentially sampling stars at lower 
heights seems to also have been slightly reduced as the bin 
size is small enough for the dispersion to be calculated using 
only local stars. A simple linear fit to the (z,a z ) plot gives 
a gradient of da z /dz = 7.87kms _1 kpc -1 . 



6 CONCLUSIONS 

We have d rawn a sample of 412 stars from the distribution 
function of lBinnevI (|2012aT ). The sample was chosen to repli- 



Figure 7. iGirard et al.l method: average mean W velocity and 
W velocity dispersion, cr z , against height above Galactic plane for 
the 100 samples of 412 stars. The points give the a verage values 
calcul ated in each bin using a method inspired bv IGirard et al.l 
(|2006h outlined in Section [5] The error bars show the standard 
deviation of the calculated values in each bin. The red dashed 
lines show the velocity moments of the thick disc calculated using 
the underlying distribution function. The dispersion is slightly 
underestimated due to the neglected thin disc contamination and 
preferentially sampling stars from lower Galactic heights, but in 
general the recovery of the truth is more successful than using 
the MB method. In particular the calculated gradient in a z is 
approximately a factor of three larger. 



cate the sample presented by MB. We performed the same 
procedure as MB to extract the mean vertical velocity and 
vertical velocity dispersion of the thick disc as a function of 
Galactic height and compared it to the known moments of 
the vertical velocity of the thick disc of the underlying dis- 
tribution. We find that the variation of the dispersion with 
z is far noisier than that found by MB. The majority of the 
error arises from the Poisson noise inherent in the limited 
sample size. We show that the observational errors in the ve- 
locities, which MB claim as the total error, cannot account 
for the deviation. A larger source of error arises from dis- 
tance errors moving stars from bin to bin and we show that 
this is of a similar order to the Poisson noise. 

A large number of samples reveals that the method sys- 
tematically overestimates the dispersion at low z and under- 
estimates it at high z. The two effects combined lead to a 
flatter curve of dispersion against z. We recalculated the 
mean velocity and velocity dispersion using a method in- 
spired by G06 which reveals the causes of these two effects: 
at low z, only fitting the wings of the distribution to re- 
move thin disc contaminants makes the distribution appear 
broader, and at high z, large bin sizes preferentially sample 
stars at lower Galactic height which in general have a lower 
velocity. The G06 method produces a much better fit to the 
expected velocity dispersion and there is approximately a 
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factor of three discrepancy in the gradient of the dispersion 
as a function of z between the results of the MB and G06 
methods. This discrepancy is not a result of more precise 
measurements or analysis but purely a result of systematics 
in the data analysis introduced by MB. We have not touched 
upon the values given by MB for the U and V dispersions 
but similar effects are expected to occur. 

The results presented here should serve as a useful 
demonstration of the expected errors and potential biases 
which arise when using a method similar to the MB method. 
We have demonstrated the need to understand the errors 
and systematics of methods which are applied to observa- 
tional data and that pseudo-samples from realistic Galaxy 
distribution functions are a useful tool in this respect. The 
effect that the biases and errors demonstr ated in this paper 
have on the dark matter mass estimates (M oni Bidin et al.l 
l2012bl : iBovv &; Tremainell2012h is beyond the scope of this 
paper. However, the results of this paper are relevant to both 
these mass estimate determinations and more generally to 
the understanding of Galactic disc kinematics. 
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